Los datos usados fueron extraídos de una base de datos históricos, representan el comportamiento de un separador bifásico (Liquido / Gas). Representan las siguientes variables:
Los datos fueron muestreados cada 2 segundos durante 2 días.
Observemos el comportamiento en una escala temporal menor:
La gráfica contiene 2500 muestras de 2 segundos, es decir 5K segundos lo cual corresponde a un intervalo de 1 hora y 25 minutos aproximadamente aproximadamente. Se observa que la señal esta sometida a ruido, pero se deja ver una comportamiento tendencial y posiblemente uno cíclico.
Se observa el comportamiento de las señal en el dominio de la frecuencia para determinar la existencia de patrones ciclicos.
Tenemos las frecuencias concentradas en la parte baja del espectro, lo cual se verifica en las siguientes gráficas:
Se observa un segundo pico en la frecuencia 0.0047 aproximadamente es decir, un periodo de aproximadamente 213 muestras, 416 segundos o 6.9 minutos.
La auto-correlación es una medida de que tanto se parece una señal a si misma cuando es desplazada en el tiempo, permite identificar patrones cíclicos con facilidad:
Se hace evidente el comportamiento cíclico de la señal, para determinar el lag efectivo para el que se tiene la periodicidad:
[1] "Lag de periodicidad: 209"
El lag 209 equivale a 418 segundos o 6.97 minutos, para efectos prácticos vamos a tomar 7 minutos.
Para poder realizar un análisis sobre la varianza real asociada a la señal de flujo, primero se debe asegurar que se segrega la tendencia y la componente cíclica, para lo cual hacemos dos diferenciaciones, la señal residual es el objeto de estudio:
[1] "Media del ruido de diferenciación: -0.0303740329780939"
[1] "Desviacion estándar de ruido de diferenciacion: 724.984381535395"
Se estudia ahora como cambia la incertidumbre de la medición, para esto se divide el dataset en paquetes de 500 muestras, para cada una de los cuales se determina la incertidumbre asociada al mejor modelo arima de cada paquete:
En el gráfico se ve la evolución de la incertidumbre de medición para los diferentes intervalos, la banda comprendida entre las lineas rojas representa el intervalo de confianza del 95% alrededor de la media de la incertidumbre de medición, los porcentajes se evalúan en comparación con el valor promedio del flujo. Nótese que durante los primeros momentos la incertidumbre de medición se encontraba muy alejada del intervalo de confianza, lo que podría indicar alguna anomalía en el funcionamiento del instrumento durante esos primeros instantes.
A continuación se reporta el error medio de medición y los limites del intervalo de confianza expresados en unidades de ingeniería:
[1] "Error de medición medio: 673.6 [BPD]"
[1] "Limites indicadores de anomalia: [200.94, 1146.26] [BPD]"
Ahora que conocemos el comportamiento en frecuencia de nuestra señal, se selecciona un filtro que elimine el ruido de la señal original observemos el sistema filtrado:
Como es el ruido si elimino la señal filtrada?
[1] "Media del ruido filtrado: 0.000547632095091429"
[1] "Desviacion estándar de ruido filtrado: 705.647413586267"
Con un filtro que hace un promedio ponderado sobre 10 muestras de la señal original, se encuentra un ruido de características similar a lo encontrado con la diferenciación (sección Análisis de incertidumbre de la medición). Este filtro sería recomendable para aplicar en la medición de flujo.
Supongamos que se desea estudiar el nivel del separador, por lo que se define este como la variable de salida de este análisis, observemos primero la correlación entre las variables:
PIT704B3 PIT704B2 PIT704B1 LIT704B1 LCV704B1 FIT704B1 LIT704B2
[1,] -0.402521 -0.4254408 -0.394431 0.8358122 0.4333321 0.2237703 1
FCV704B3 LIT701A PCV704B1
[1,] 0.3577001 0.06468231 0.1833118
Se observa que muchas variables tienen altas correlaciones con el LIT704B2, que es la variable de nivel de interés, obsérvese, por ejemplo, que la correlación con el LIT704B1 es de mas del 85%, lo cual se explica por el hecho de que el LIT704B2 es una redundancia del LIT704B1 (de hecho, es curioso que la correlación no sea mayor). Las correlaciones ademas indican que la menor correlación se presenta con LIT701A, que representa el nivel en el tanque que esta aguas abajo del separador.
Teniendo en cuenta la forma en que funciona el separador bifásico, también llama la atención que la correlación entre el nivel y el porcentaje de apertura de la válvula de salida no sea mayor, verifiquemos la cros-correlación entre estas dos señales:
La cros-correlación muestra la correlación de una señal con otra señal en distintos corrimientos en tiempo, obsérvese como hay un patrón cíclico para el que se tiene que la cros-correlación alcanza nuevos máximos de forma periódica, determinemos cual es es periodo de tiempo:
[1] "Numero de muestras por periodo: 217"
[1] "Tiempo muestras por periodo: 434 [seg]"
El periodo de tiempo para cada repetición es de aproximadamente 7.23 minutos, comparemos el comportamiento de la válvula con el comportamiento del nivel:
Se tiene que ambas señales se mueven de forma muy parecida, pero están desfasadas, la señal de la válvula (en rojo) toma acción mucho después de que ocurra un cambio en el nivel, lo anterior indica que hay espacio para la mejora en la sintonización del lazo de control o quizá hay un problema en la válvula que hace que responda tardíamente a los cambios en el nivel. De la gráfica de cros-correlación podemos determinar cual es el retraso en la respuesta de la válvula, en el caso ideal, las señales observadas deberían tener la mayor cros-correlación para el corrimiento cero, pero como están desfasadas la mayor cros-correlación se presentará en el corrimiento de desfase:
[1] "Numero de muestras de desfase: 25"
[1] "Tiempo de retardo de respuesta de la válvula: 50 [seg]"
Dentro del análisis de un sistema se procura eliminar las señales redundantes, ya que estas inducen a errores en la determinación de la verdadera calidad del modelo desarrollado, la redundancia se observa calculando el coeficiente de inflación de varianza:
PIT704B3 PIT704B2 PIT704B1 LIT704B1 LCV704B1 FIT704B1 FCV704B3 LIT701A
2.861677 74.696278 74.156571 1.591187 2.813205 2.169959 3.004261 1.255956
PCV704B1
1.482531
La inflación de varianza se considera normal para valores por debajo de 5, obsérvese que los PIT704B1/B2 superan los valores normales de inflación de varianza, lo cual es consistente con el hecho de ser instrumentos redundante, LIT704B1 es también un instrumento redundante, se deben eliminar dichas redundancias para determinar el desempeño real del modelo:
PIT704B3 PIT704B2 LCV704B1 FIT704B1 FCV704B3 LIT701A PCV704B1
2.744941 1.109100 2.462284 2.155562 2.896113 1.214823 1.463190
Observemos como son las relaciones entre el nivel en el separador y las demas variables gráficamente:
Es difícil encontrar alguna relación clara, sin embargo, PIT704B2, LCV704B1 y FIT704B1 parecen tener algún comportamiento sistemático con relación al LIT704B2, desarrollemos un modelo usando regresión lineal para determinar mas claramente las relaciones entre las variables:
Call:
lm(formula = LIT704B2 ~ ., data = df_w1)
Residuals:
Min 1Q Median 3Q Max
-1.81604 -0.35075 -0.00244 0.33806 2.37416
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.054e+05 3.078e+03 -34.231 < 2e-16 ***
date 6.534e-05 1.901e-06 34.378 < 2e-16 ***
PIT704B3 -1.263e-01 6.034e-03 -20.929 < 2e-16 ***
PIT704B2 -1.185e+00 3.519e-02 -33.672 < 2e-16 ***
LCV704B1 3.201e-01 3.590e-03 89.170 < 2e-16 ***
FIT704B1 -5.802e-06 2.096e-06 -2.768 0.005640 **
FCV704B3 1.591e+00 4.571e-01 3.481 0.000501 ***
LIT701A -1.440e+00 4.350e-02 -33.099 < 2e-16 ***
PCV704B1 9.189e-02 2.426e-03 37.880 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5348 on 14992 degrees of freedom
Multiple R-squared: 0.56, Adjusted R-squared: 0.5597
F-statistic: 2385 on 8 and 14992 DF, p-value: < 2.2e-16
El modelo explica aproximadamente el 56% de la variabilidad de la señal de nivel, observemos como son los residuos de la regresión (es decir, aquella parte no puede explicar el modelo)
Los errores residuales parecen mostrar un comportamiento aleatorio normal, casi siempre con incertidumbres menores a +/- 1 pulgada, sin embargo, por la naturaleza misma del proceso, se sabe que hay grandes posibilidades de que el ruido residual este correlacionado y sea una función del tiempo, esto resulta en que la relevancia real de las variables se ve enmascarada, dando cabida a darle importancia a variables que realmente no la tienen.
Adicionalmente, debido a la gran cantidad de muestras disponibles, las pruebas estadísticas de los datos correlacionados suelen arrojar p-values demasiado pequeños. Analicemos el sistema tomando solo 4000 muestras:
Call:
lm(formula = LIT704B2 ~ . - date, data = df_w1)
Residuals:
Min 1Q Median 3Q Max
-1.44442 -0.28512 0.03313 0.27903 1.36918
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.964e+02 6.533e+01 6.067 1.42e-09 ***
PIT704B3 2.024e-02 9.053e-03 2.236 0.0254 *
PIT704B2 -1.048e+00 5.261e-02 -19.918 < 2e-16 ***
LCV704B1 3.952e-01 6.481e-03 60.989 < 2e-16 ***
FIT704B1 -9.557e-06 2.817e-06 -3.392 0.0007 ***
FCV704B3 -3.335e+00 6.524e-01 -5.111 3.35e-07 ***
LIT701A -3.193e+00 8.943e-02 -35.704 < 2e-16 ***
PCV704B1 9.083e-02 3.974e-03 22.856 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4014 on 3993 degrees of freedom
Multiple R-squared: 0.6827, Adjusted R-squared: 0.6821
F-statistic: 1227 on 7 and 3993 DF, p-value: < 2.2e-16
Box-Pierce test
data: fit1$residuals
X-squared = 26986, df = 8.2943, p-value < 2.2e-16
En este intervalo de tiempo es mucho mas claro que el ruido esta correlacionado (el valor del ruido ) La prueba de Box-Pierce sirve para determinar si las auto-correlaciones de los residuos son cero (es decir no hay), si el p-value es menor a 0.05 se rechaza la hipótesis (de que no hay auto-correlación). Dado que el test dio un p-value mucho menor a 0.05 existe fuerte evidencia de que el ruido tiene auto-correlaciones.
La auto correlación parcial sugiere un comportamiento autoregresivo de grado 2 o 3, Observemos ahora el ACF para mas retardos:
Es evidente que se tiene un remanente de comportamiento cíclico, sobretodo teniendo en cuenta la presencia de ciclos repetitivo cada 7 minutos, revisemos la densidad espectral:
Calculemos los máximos para la densidad espectral:
# A tibble: 6 x 4
freq spec per_seg per_min
<dbl> <dbl> <dbl> <dbl>
1 0.00469 47.4 426. 7.11
2 0.00444 28.4 450 7.5
3 0.00494 24.5 405 6.75
4 0.00519 16.8 386. 6.43
5 0.00321 16.3 623. 10.4
6 0.00148 15.7 1350 22.5
se observa de nuevo que la mayor potencia se concentra en un periodo cíclico de aproximadamente 7 minutos, existen otras componentes que pueden ser significativas en 10 y 22 minutos, sin embargo no las consideraremos para el ejercicio.
Consideremos diferenciar los residuales y comparar el ruido obtenido:
El ruido final parece ser mucho mas estacionario, se ha eliminado la tendencia remanente y el patrón cíclico (obsérvese que en la figura anterior se muestran las dos fases, al residuo se le extrae la parte cíclica y se obtiene la curva en rojo, luego a este se le extrae la tendencia y se obtiene la curva verde)
Finalmente al resultado se le llama proceso estacionario y se gráfica a continuación:
Grafiquemos nuevamente la auto-correlación y la auto-correlación parcial:
Parece que podemos tener un grado 2 AR y un grado 2 AM, para determinar que modelo se ajusta mejor simulemos todas las posibilidades y escojamos la mejor:
ar_opt ma_opt
2 0
Se obtiene que el modelo que mejor resultado ofrece es AR grado 2, después de haber obtenido todos estos resultados, podemos construir el modelo para el sistema completo (el resultado mostrado es producto de eliminar las variables que no son representativas estadisticamente hablando):
Estimate SE t.value p.value
ar1 0.1665 0.0259 6.4255 0e+00
ar2 -0.0982 0.0246 -3.9870 1e-04
LCV704B1 0.0398 0.0056 7.1477 0e+00
LIT701A 0.3174 0.0914 3.4732 5e-04
PCV704B1 -0.0652 0.0072 -9.0175 0e+00
Los residuales resultantes de la regresión se muestran en las gráficas a continuación:
Los residuales del modelo resultante muestran aun comportamientos cíclicos, esta vez con un periodo de 4.4 minutos aproximadamente (ver gráfica de auto-correlación, el periodo es de 132 muestras).
Del análisis desarrollado se pueden resumir varios hallazgos relacionados con el sistema estudiado: